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Abstract 

The  acoustic  wavefield  reflection  response  of  a  plane  ware  on  a 
cylindrical  surface  is  calculated  from  a  specialization  of  the  Eirchhoff 
integral.  Conputational  advantages  are  obtained  by  assuming  that  structural 
changes  occur  only  along  the  direction  of  data  collection.  Off-line 
geologic  invariance  permits  an  integral  to  be  replaced  by  a  short 
convolution  operator.  The  restriction  to  single  interface  modeling  permits 
implementation  on  personal  computers.  Also,  the  single  layer  algorithm 
provides  the  framework  for  a  multi-layer  code. 
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Introduction. 


An  efficient  program  has  been  developed  to  model  the  2-dimensional 

Kirchhoff  acoustic  vavefield  reflection  response  of  a  plane  wave  on  a 
cylindrical  surface.  Intuitively,  a  source  function  is  placed  on  each 
differential  surface  element,  and  an  exploding  reflector  simulation  is 

performed  by  initiating  all  of  the  sources  concurrently.  Ihe  method 
discussed  here  is  based  on  a  technique  by  John  R.  Berryhill  (1979)  in  which 
cylindrical  symmetry  allows  for  considerable  execution  time  reduction  in 
computing  the  Rayleigh  II  integral.  In  brief,  assuming  invariance  of 

geologic  structure  in  the  off-line  direction  permits  a  short  convolution 
operator  to  replace  one  of  the  two  requisite  surface  integrals. 

There  is  an  advantage  in  having  a  program  specialized  to  single-layer 
modeling  in  that  the  normal-incidence  wavefield  effects  of  a  plane  wave  on 
buried  anticlines,  synclines,  or  arbitrary  structures  may  be  simulated 
without  the  computational  overhead  associated  with  performing  single-layer 
responses  via  a  multi-layer  program.  The  small  memory  requirements  and 

relatively  low  amount  of  input/output  make  the  program  particularly  friendly 
to  implementation  on  a  personal  computer  (80  traces  in  10  minutes  on  an  IBM 
PC/XT  with  an  8087  co-processor).  Also,  the  single  layer  modeling  routine 
is,  hopefully,  a  large  step  toward  developing  a  multi-layer  code. 


A  general  description  of  the  Berryhill  technique  as  applied  to  wave 
equation  detuning  is  given.  For  completeness,  a  derivation  of  the  algorithm 
is  included.  Example  simulations  are  then  provided  which  demonstrate  the 
wavefield  effects  inherent  in  the  method. 


-  2  - 


Wave  Equation  D>twii| 


The  primary  application  of  the  extrapolation  method  employed  by 
Berryhill  (1979)  is  in  wave  equation  datnming.  This  procedure  redefines  the 
surface  on  which  the  data  are  taken,  for  the  purpose  of  minimizing 
distortive  surface  effects.  A  common  example  inclndes  downward  continuing 
marine  data  to  an  irregular  ocean  bottom,  and  then  upward  continuing  the 
data  to  a  flat  plane,  using  the  velocity  from  under  the  ocean  bottom.  This 
removes  the  effect  that  an  irregular  surface  has  on  deeper  reflectors. 

In  a  modeling  sense,  a  source  is  placed  on  a  cylindrical  surface  and 
the  resulting  wavefield  is  recorded  at  another  surface.  Note  that  only  an 
inter-datum  velocity  is  required.  The  technique  is  therefore  based  on  the 
'exploding  reflector"  model  of  the  subsurface  using  a  velocity  equal  to  half 
the  true  velocity  in  order  to  produce  correct  arrival  times. 

In  a  multi-layer  setting,  the  initial  exploding  reflector  corresponds 
to  the  deepest  horizon,  with  an  identical  source  placed  on  each  differential 
surface  element.  The  resulting  wavefield  is  then  recorded  at  the  next 
reflecting  surface.  A  time  section  displaying  the  wavefield  as  seen  by  the 
upper  surface  is  thus  produced.  To  each  trace  in  this  time  section,  the 
original  source  function,  weighted  by  a  reflection  coefficient,  is  added. 
When  this  datum  explodes,  it  extrapolates  the  effects  of  the  lower 
reflector,  true  to  Snell's  law,  as  well  as  extrapolating  its  own  exploding 
datum  upwards.  The  same  procedure  is  applied  to  each  reflector,  as  the 


Algorithm  Derivation 


The  Berryhill  algorithm  is  an  application  of  the  Rayleigh  II 
specialization  of  the  Kirchhoff  integral.  The  following  ontline  of  the 
Rayleigh  II  integral  derivation  is  from  Berkhout  (1982).  Fnrther  details 
can  be  found  in  texts  by  Born  and  Volf  (1983),  and  Lipson  and  Lipson  (1981). 

Define  P  as  the  pressure  inside  a  region,  bounded  by  a  surface  S,  due 
to  a  compressional  source  exterior  to  S.  The  free-space  Green's  function, 
G,  represents  a  pressure  due  to  an  interior  source.  Therefore 


and 


+  £*)  P  =  0  . 

+  K2)  G  =  -  4jt8(x-x  )8(y-y  )8(z-z  ) 

ft  8  ft 


Applying  Green's  theorem  results  in 


(( PVG  -  GVP)  •  n  ds  =  -  4n  P(x  ,y  ,z  )  =  -  4n  P 
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Noting  that 
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where  pQ  is  the  static  density  and  VQ  is  the  normal  component  of  particle 
velocity  leads  to  the  Kirchhoff  integral. 
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ds 


Vhere  A  is  now  defined  as  the  inward  pointing  normal,  and 


cosp  =  - 

as 

The  Kirchhoff  integral  states  that  the  pressure  at  a  point  A  is  a  surface 
extrapolation  of  monopole  normal  velocity  components,  and  a  dipole  pressure 
field. 


The  surface  is  chosen  to  be  a  closed  hemisphere,  such  that  the  distance 
between  the  point  A  and  the  planar  portion  of  the  hemisphere  is  much  less 
than  the  distance  between  the  point  A  and  the  acnrate  portion.  The 
integration  surface  is  therefore  the  planar  base  of  a  very  large  hemisphere. 
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An  alternate  (negative  image)  Green's  function  eliminates  the  normal 


velocity  component  term,  providing  the  Rayleigh  II  integral. 
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It  is  interesting  to  note  that  applying 


dP  3V 
_____  __  _ n 

In  Po  at 

to  the  Rayleigh  I  integral  produces  the  Rayleigh  II  integral,  vith  P 
replaced  by  Vn.  Thus,  the  Rayleigh  II  integral  can  be  aeen  as  a  means  of 
forward  (or  reverse)  extrapolation  of  either  a  dipole  pressure  field,  or  the 
normal  components  of  particle  velocity. 

Berryhill  (1979)  starts  with  the  time  domain  representation  of  the 
Rayleigh  II  integral,  and  then  assumes  that  there  is  no  variance  of  the 
extrapolated  surface  in  the  y-direction.  Under  this  assumption,  and  in 
Berryhill 's  notation. 
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where  the  geometry  is  defined  in  Fignre  1. 


The  y-integration  is  replaced  with  a  convolntion,  leading  to 


00 

_  1  I  c  du  ( x,  o,  t )  *  „  ,  ^  ,  v  d  tan  r 
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»  I  cr^  dt  x  dt 

where  H(t)  is  the  unit  step  function,  and 
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In  the  digital  implementation  of  the  algorithm,  the  derivative  on  the 
source,  u(x,  o,  t),  is  moved  to  the  tangent  operator. 


This  new  tangent  operator,  containing  the  off-line  integral 
contributions,  goes  to  zero  rapidly.  Typically,  the  first  9  points  of  the 
operator  do  nicely,  for  computer  results,  as  seen  in  Fignre  2.  Therefore, 
the  response  at  a  point  due  to  a  wavefield  defined  on  a  surface  is  a  sum  of 
line  responses  on  the  surface.  A  single  line  response  is  illustrated  in 
Figure  3.  Each  line  response  is  thus  obtained  by  convolving  the  surface 
source  function  with  a  short  off-line  operator. 


A  final  assumption  is  made  in  defining  a  curved  surface  as  a  sequence 
of  properly  weighted  planar  segments.  The  diffraction  response  due  to 
segment  corners  is  negligible  if  the  associated  dihedral  angles  are  less 
than  a  few  degrees. 
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CONVOLUTION  OPERATOR  FOR  Y-I NTEGRATION 


FIGURE  2 


FIGURE  1 

Following  Berryhill  [19791 
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kept  constant,  with  the  height  of  the  feature  increasing.  The  diffraction 
energy  in  the  corresponding  tine  sections  illustrates  the  fidelity  of  the 
method  of  producing  vavefield  responses. 
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Examples 

All  of  the  examples  in  Figure*  10-17  were  generated  using  the  flat- 
topped  averaging  operator.  The  ont  of  plane  operator  length*  n*ed  were 
extravagantly  long,  and  computationally  wasteful.  A  maxima  of  10  points 
works  lust  as  well/  (do  as  I  say,  and  not  as  I  dol).  Receivers  are  on  a 
flat  datum,  bnt  need  not  necessarily  be  so.  A  cubic  spline  is  fit  to  both 
the  receiver  and  reflector  surfaces.  The  two  lines  drawn  to  the  exploding 
surface  from  each  receiver  depicts  the  boundaries  of  the  respective 
integration  zones. 

A  search  is  conducted  for  spectral  points  associated  with  each 
receiver.  The  search  zone  is  symmetric  about  the  receiver,  and  8  Fresnel 
zones  wide.  The  location  of  an  extreme  right  and  an  extreme  left  spectral 
point  delineate  the  integration  zone,  the  outer  borders  of  which  are  one 
Fresnel  zone  beyond  each  spectral  point.  If  no  spectral  point  is  found,  the 
integration  is  performed  in  a  zone  directly  under  the  receiver. 

Figure  10  illustrates  the  geometry  of  a  syncline  model.  Receivers  are 
positioned  along  the  top  of  the  diagram,  with  their  respective  integration 
zones  delineated.  The  resultant  time  section  in  Figure  11  depicts  the  "bow 
tie”  effect  associated  with  synclinal  structures.  Note  that  the  waveform 
appearing  in  the  center  of  the  bow  tie  is  the  derivative  of  the  Ricker 
source  function.  This  is  due  to  the  buried  focus. 

Fugure  12  through  17  model  the  diffraction  effects  of  an  anticline. 
The  distance  between  the  receiver  datum  and  the  base  of  each  structure  is 
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increment  is 


s  function  of  a  maximum  frequency. 
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frequency  domain  aliasing  is  displayed  via  a  standard  scale  and  a  zoom  plot. 

Tlie  averaging  operator  permits  a  constant  sample  rate  along  the 
integration  surface  by  filtering  out  the  high  frequency  components  of  the 
late  arrivals.  This  is  accomplished  by  convolving  each  segment  response 
with  a  flat-topped  averaging  operator.  The  length  of  the  operator  being 
proportional  to  the  difference  in  arrival  times  of  adjacent  segments. 
Figure  7  displays  the  individual  segment  contributions  after  application  of 
the  averaging  operator.  Note  that  the  frequency  content  of  the  late  arrivals 
has  been  reduced.  Summation  of  the  line  sources  results  in  a  much  cleaner 
vavelet,  as  displayed  in  Figure  8. 

The  source  function  in  all  examples  is  a  causal  Ricker  vavelet.  As 
depicted  in  Figure  9,  it  is  symmetric  in  the  time  domain,  consisting  of  a 
small  trough,  a  large  peak,  and  a  small  trough.  In  the  frequency  domain, 
the  vavelet  is  approximately  Gaussian  shaped,  vith  "important”  information 
coming  from  frequencies  of  up  to  2  *  peak  frequency.  The  Nyquist  criteria 
vas  therefore  based  on  2  *  peak  frequency. 

The  area  over  vhich  to  perform  the  integration  is  crucial.  In  terms  of 
a  Fresnel  zone,  the  integration  zone  should  be  as  vide  in  space  as  the 
source  is  long  in  time.  The  "important”  lov  frequency  determines  the  large 
vavelength  for  calculating  the  Fresnel  zone  vidth.  Typically,  a  value  of 
1/2  *  peak  frequency  suffices.  (It  should  be  noted  that  taking  advantage  of 
cylindrical  symmetry  has  dictated  a  rectangular  "Fresnel  zone”,  instead  of 
the  usual  circle.)  Therefore,  the  integration  zone  vidth  is  determined  by  a 
minimum  frequency,  and  for  a  given  sub-surface  integration  area,  the  surface 


Iaplmititioa 


A  numerical  integration  is  carried  ont  along  the  exploding  datum 
employing  a  constant  X  increment.  If  the  distance  between  segments  is  too 
large  (relative  to  "important*  short  wavelengths)  aliasing  occurs, 
manifesting  itself  in  the  form  of  very  *ringy*  late  arrivals.  A 
convolutional  averaging  operator  may  be  applied  if  desired.  and  if 
necessary.  If  sampling  is  done  properly.  the  averaging  operator 
automatically  reduces  to  a  delta  function,  and  thus  has  no  effect. 

For  each  receiver  location  a  Fresnel  zone  is  calculated,  proportional 
to  the  longest  important  wavelength  in  the  source  fnnotion.  The  numerical 
spatial  sampling  is  then  determined  such  that  the  high  frequency  content  of 
the  late  segment  contributions  are  not  aliased. 

As  detailed  previously,  the  recorded  response  from  an  exploding  plane 
is  obtained  by  summing  responses  of  exploding  lines.  Figure  4  illustrates 
an  exploding  planar  reflector  at  the  z  =  5000  foot  level  and  a  receiver  at 
the  z  =  10,000  foot  level.  Each  line  drawn  from  the  plane  to  the  receiver 
represents  a  single  line  source  contribution  in  the  plane.  The  individual 
contributions  before  summation  are  plotted  in  Figure  5.  As  long  as  the 
spatial  sampling  is  fine  enough,  when  the  contributions  are  added  together 
they  should  produce  a  wavelet  resembling  the  Ricker  source  wavelet.  In 
Figure  6  it  is  seen  that  the  sampling  is  adequate  for  that  portion  of  the 
surface  directly  under  the  receiver  location.  However,  as  the  distance 
increases  to  each  segment,  so  does  the  travel  time,  and  eventually 
contributions  are  arriving  too  far  apart,  and  aliasing  occurs.  In  the 
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